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We present a theoretical study on the orbital magnetism in multilayer graphenes within the 
effective mass approximation. The Hamiltonian and thus susceptibility can be decomposed into 
contributions from sub-systems equivalent to monolayer or bilayer graphene. The monolayer- type 
subband exists only in odd layers and exhibits a delta-function susceptibility at ef — 0. The bilayer- 
type subband appearing in every layer number gives a singular structure in the vicinity of ef — 
due to the trigonal warping as well as a logarithmic tail away from ef = 0. The integral of the 
susceptibility over energy is approximately given only by the layer number. 



I. INTRODUCTION 

Recently unconventional electronic properties of mono- 
crystalline graphenes attracts much attention motivated 
by experimental fabrication)^^ although they were al- 
ready the subject of theoretically study prior to the 
fabrication4^^^iii^iH Multilayer films which con- 
tain more than two layers can also be synthesized, and 
various phenomena depending on the layer number have 
been reported^ ' 14 ' 15 In this paper we present a the- 
oretical study on the orbital magnetism in multilayer 
graphenes. 

The electronic structure of the monolayer graphene 
is quite different from conventional metals, because 
the conduction and valence bands touch at K and K 1 
points in the Brillouin zone, around which the disper- 
sion becomes linear like a relativistic particle. In multi- 
layer graphenes, the interlayer coupling makes a com- 
plex structure around the band touching. The elec- 
tronic properties of graphene bilayer were theoretically 
studied for the band structur e 16 ! 17 and the transport 
properties. 18 ' 19 ! 20 For few-layered graphenes of more than 
two stacks, the electronic structure is investigated theo- 
retically in a k • p approximation^ a density functional 
calculation^ 2 , and a tight-binding modeh 23 ! 24 On the ex- 
perimental side, the band structures of graphenes from 
one to four layers were recently measured using angle- 
resolved photoemission spectroscopy^ 

The orbital magnetism in graphene-based systems was 
first studied for a monolayer as a simple model to ex- 
plain the large diamagnetism of graphite X It was found 
that the susceptibility becomes highly diamagnetic at 
e = (band touching point) even though the density 
of states vanishes there. The calculation was extended 
to graphit o 25 ' 26 and to few-layered graphenes as a model 
of graphite intercalation compounds! 27 ' 28 ' 29 The Fermi 
surface of the graphite is known to be trigonally warped 
around the band touching point^ 5 - and the effect of the 
warping on magnetization was discussed within the per- 
turbational approach, 26 Recently, the disorder effects on 
the magnetic oscillatio n 30 ' 32 and on the susceptibilit y 31 ' 32 
were studied for the monolayer graphene. 



Here we present a systematic study on the orbital 
magnetism for multilayer graphenes with arbitrary layer 
numbers in the effective mass approximation. We show 
that the Hamiltonian of a multilayer graphene can be de- 
composed into those equivalent to monolayer or bilayer, 
which allows us to study the dependence of the suscep- 
tibility on layer numbers. We take the trigonal warping 
effect into the calculation and show that the fine struc- 
ture around zero energy gives rise to singular magnetic 
properties. 

We introduce the model Hamiltonian and its decom- 
position into subsystems in Sec. HH and present the cal- 
culation of the magnetization in Sec. IIIIl The discussion 
and summary are given in Sec. IIVI 



II. FORMULATION 

We consider a multilayer graphene composed of N lay- 
ers of a carbon hexagonal network, which are arranged 
in the AB (Bernal) stacking, as shown in Fig. [T] A unit 
cell contains Aj and Bj atoms on the layer j = 1, • • • , N. 
For convenience we divide carbon atoms into two groups 
as 

Group I: B x , A 2 , B 3 , ■ ■ ■ (I) 
Group II: A U B 2 ,A 3 ,--- (2) 

The atoms of group I are arranged along vertical columns 
normal to the layer plane, while those in group II are 
above or below the center of hexagons in the neighboring 
layers. The lattice constant within a layer is given by 
a = 0.246 nm and the distance between adjacent layers 
c /2 = 0.334 nm. 

The system can be described by a kp Hamiltonian 
closely related to a three-dimensional (3D) graphite 
modeh 25 ' 33 ' 34 ' 35 The low energy spectrum is given by the 
states in the vicinity of K and K' points in the Bril- 
louin zone. Let \Aj) and \Bj) be the Bloch functions 
at the K point, corresponding to the A and B sublat- 
tices, respectively, of layer j. For monolayer graphene, 
the Hamiltonian around K point for the basis \B{) 
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with 70 being the nearest-neighbor coupling in a single 
layer. We cite the experimental estimation 70 w 3.16 
eV^ 

For the inter-layer coupling, we include parameters 
71, and 73, following the Hamiltonian previously derived 
for a bilayer graphene i 16 i 18 Here 71 represents the cou- 
pling between vertically neighboring atoms in group I 
(A 2 k «-> B2k±i), and 73 between group II atoms on neigh- 
boring layers (B 2 k <-» A 2 k±i), which are estimated to 
71 « 0.39 eVi^ and 73 » 0.315 eV^i If we look at the 
interaction between layers 1 and 2, the matrix element 
(A 2 \H\Bi), corresponding to the vertical bond, becomes 
71 not accompanying in-plane Bloch number. The matrix 
element (B 2 \H\Ai) is written as r y'k + with 
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-A73, 



(5) 



similar to the intra-layer term (Ai\H.\Bi) , as the in-plane 
vector components from A\ to B 2 are identical to those 
from B\ to A\, 

Accordingly, if the basis is taken as |Ai),|£>i); 
\A 2 ), I-B2); • • • ; |Aiv)) I-Bjv), the Hamiltonian for the mul- 
tilayer graphene around the K point becomes 
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with 



V = 



7 'fc 4 
7i 



The effective Hamiltonian for K' is obtained by exchang- 
ing k + and fc_ and replacing 71 with —71. The deriva- 
tion of the effective Hamiltonian based on a tight-binding 
model is presented in Appendix lAl 

We show in the following that the Hamiltonian matrix 
(|6|) can be block-diagonalized into smaller matrices by 
choosing an appropriate basis independent of k. First, 
we arrange the basis in the order of group I and then 
group II, i.e., \B X ), \A 2 ), \B 3 ), • • • ; \A{), \B 2 ), \A 3 ), 
Then, Eq. ([6]) becomes 



\H\ 2 H 22 
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FIG. 1: Atomic structure of multilayer graphene with AB 
(Bernal) stacking. 



with Hij being N x N matrices defined as 
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where the upper and lower signs correspond to odd and 
even TV, respectively. 

If we set k = (the K point), Tti 2 and Tt 22 van- 
ish. Remaining Tin is equivalent to the Hamiltonian of 
a one-dimensional tight-binding chain with the nearest- 
neighbor coupling 7! , giving a set of eigenenergies 



= 7lAjV,-, 

= 2 sin 



2(N + iy 



with 



m= —(N - 1), -(iV-3), 



N — 1. 
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Here, m is an odd integer when the layer number N is 
even, while m is even when N is odd, and therefore m = 
is allowed only for odd N. 

The corresponding wave function is explicitly written 



as 



V'mO') = 



N + l 



■ sin 



(-m + N + l)n .' 
2(N + 1) 3 



(14) 
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where tp m (j) represents the amplitudes at \B\), \A 2 ), 
\B 3 ), ■ ■ ■ and satisfies 



^2^7n(j)lp7n'(j) = S n 



(15) 



We have a relation between the wave functions xp m and 

■0-m as 



3+1 



(16) 



Now we construct the basis by assigning ip m (j) to the 
atoms of group I and II as 



(17) 



= V>m(l)|£l) + 1pm(2)\A 2 ) + ^m(3)|B 3 ) + ' ■ 
= ^m(l)|Ai) + ^ m (2)|B 3 ) + ^„(3)|A 3 ) 4 



and attempt to rewrite the Hamiltonian ([6]). The ma- 
trix elements within group I come from TL\\ and become 
diagonal as is obvious from the definition, 



(18) 



Off-diagonal elements between |0m ) and are writ- 

ten from 7ii2 as 

<*£?l«ltf£?> 

N N 

3=1 3 = 1 

= ^{kx&rn.m' iky$m,—m')- (19) 

In the second equality we used relation (flB"]) and orthog- 
onality (|15[) . Lastly, the matrix elements within group II 
are obtained from H22 as 



^ hCO' + l)^m(j) + C'(j')^m(j + 1)] 



3=1 

+y<*» E [^O' + i)^(j) - r m >uwm(j + 1)] 

3=1 

). (20) 

The Hamiltonian is thus closed in the subspace 

{|0m>, l^l), |0£°>, 10-™)} for each |m|. Particularly, 
m = is special in that the subspace is spanned with 
only two bases {I^q 11 ^), |</>o^)}, while this is absent in 
even-layer graphenes. The sub-matrix is written as 



T~i m =o — 



7fc_ 

iK or 



(21) 



which is independent of 71 and 73, and equivalent to the 
Hamiltonian of the monolayer graphene. 
For to ^ 0, we rearrange the basis as 

{(iffi) + l^»/V2, + \^ m ))/V2, 

(l^?> - \<f>-l))/V2, (|^)-|^»/V2}, (22) 



where we take to > without loss of generality. We then 
obtain 



Hi, 



( 7 /c_ \ikj\ 

jk + A71 

A71 7 fc_ 

\A 7 'fc_ 7 /c+ 



(23) 



with A = \N,m- This is equivalent to the Hamiltonian 
of a bilayer graphene except that 71 and 7' (oc 73) are 
multiplied by A. 

Thus the Hamiltonian of odd-layered graphene is com- 
posed of one monolayer- type and (N — l)/2 bilayer- type 
subbands while that of even-layered graphene is com- 
posed of N/2 bilayers but no monolayer. The similar idea 
was previously proposed for trilayer graphene without 73, 
where it was shown that the energy spectrum becomes a 
superposition of that for a monolayer and for a bilayer. 21 
Here we have extended this argument to decomposition 
of the Hamiltonian matrix, and to systems with arbi- 
trary number of layers including the trigonal warping. 
We also note that k-independence of the basis becomes 
important in the following sections, since this enables us 
to write the magnetization as a sum over contributions 
from sub-Hamiltonians, which are independently calcu- 
lated. 

Many other parameters were introduced for the de- 
scription of the band structure of bulk graphite j 25 i 26 ' 35 i 38 
The parameter 74 couples group I and II atoms sit- 
ting on the neighboring layers, such as Aj <-» Aj+i or 
Bj <-> Bj + \. This parameter does not change the quali- 
tative feature of the low-energy spectrum and therefore is 
not important. 25 Parameters 72 and 75 represent vertical 
hoppings between the second-nearest neighboring layers 
for group II and I atoms, respectively. Further, 7g is an 
energy difference between the group I and II atoms due to 
difference in the chemical environment. Inclusion of these 
parameters 72, 75, and 76 causes opening up of small en- 
ergy gaps between the conduction and the valence bands. 
However, these gaps do not play important roles in the 
magnetization as will be discussed in the following. 

In 3D limit, N — > 00, the eigenstate becomes a su- 
perposition of opposite traveling waves with ±fc 2 along 
the stacking direction. The relation between the index 
to and \k z \ is obtained by comparing the eigenenergy of 
Hu, Eq. $T2j), to that of the 3D limit, 271 cos(fc z c /2), 
as 
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(-to + N + 1)tt 
2(JV+ 1) 



(24) 



The band structure of the Hamiltonian (|23|) can be 
obtained by replacing 71 by A71 and 73 by A73 in that of 
the bilayer J£ We plot in Fig. [5] the dispersion for A = 2, 
which has the maximum trigonal warping. The middle 
two subbands stick together at e = while the remaining 
two bands appear only in the energy range |e| > A71. If 
we neglect 73, the effective Hamiltonian for |e| <C A71 
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FIG. 2: (Top) Projected band structure of the sub- 
Hamiltonian Eq. (|23p with A = 2 and 73/70 = 0.1. |e| = 
£trig = O.O271 is shown as horizontal dotted lines. Right panel 
shows zoom out of the left. (Bottom) 3D plot of the lower 
second band around the band touching point. Four Fermi 
points at e = indicated by dots. 



becomes 
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) + \<t>W))/V2, {\^)-\<P (n i))/V2\, (26) 

giving a rotationally symmetric dispersion with the effec- 
tive mass 



which works for the reduced basis 
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The term proportional to 73 is responsible for the trig- 
onal warping effect, which is most remarkable around the 
band sticking point e = 0. Let us define 



A73 
7o 



(A7i)- 



(28) 



In the energy range |e| <£trig, the Fermi line splits into 
four separated pockets, one center part and three leg 
parts located trigonally, which shrink into four Fermi 
points linearly with e — > 0. We note that etrig is propor- 
tional to A 3 and thus very sensitive to A, while the en- 
ergy of the higher-band bottom, A71 , behaves linear to A. 
The maximum of A approaches 2 as the layer number in- 
creases, so that e trig becomes as large as 2(7 3 /7 ) 2 7i ~ 8 
meV. 



Figure[3]shows the band structures around the K point 
along the k x axis in the multilayer graphenes with N = 2, 
3, 4, and 5 and 73/70 = 0.1. The lists of Xn,™ are given 
as 

N = l: {A 1>0 } = {0} 

N = 2: {A 2>1 } = {1} 

A = 3: {A 3 , , A 3 , 2 } = {0,V2} 

N = 4 : {A44 , A 4 , 3 } ={(VE- l)/2, (V5 + l)/2} 

A = 5: {A 5 ,o, A 5 , 2 , A 5j4 } = {0,1,\/3}. (29) 

While we have included 70 , 71 , and 73 in our graphcne 
model, the extra parameter neglected here may make 
some changes in the electronic structure. The energy 
band of a few-layered graphene has been calculated in 
the density functional calculation 22 and the tight-binding 
modeL 23 ' 24 Those results differ from ours mainly in that 
the band centers relatively shift depending on m, and 
that a narrow gap opens where the conduction and va- 
lence bands (within a single to) touch, and where differ- 
ent bands (with different m's) cross. Gaps are attributed 
to effects of couplings such as 72, 75, and 76, which are 
mentioned above. In terms of the effective mass Hamil- 
tonian (O, those parameters appear as matrix elements 
without being multiplied by the wave number k x and 
k y , since they are associated with a hopping along the z 
axis or a diagonal element. Thus, they do not vanish at 
k = (K or K 1 ) and lift the degeneracy to open a gap. 
Apart from the gap opening, the main feature of the trig- 
onal warping is well described in the present model. It 
should also be mentioned that an energy gap is induced 
by an electric field perpendicular to the layer stacking 
direction ) 15 ! 17 ! 21 ' 42 where the electrostatic potential ap- 
pears as matrix elements independent of k x and k y as 
well. 



III. MAGNETISM OF MULTILAYER 
GRAPHENES 

For the magnetic susceptibility, we use the general ex- 
pression based on the linear response theory^ 



Im 



def(e)F(e + iO), 



(30) 



with 

F{z) 



e 2 



■1-1 2 h 2 ^ 

k 



tr ( GH X GHy GTL X GH y ) , (31) 



where g v — 2 is the valley degeneracy, g s = 2 is the spin 
degeneracy, and L is the system size. We defined here 
H x = dH/dk x , Hy = dH/dky, G{z) = (z - H)' 1 , and 
/(e) = [1 + e (£-M)/fcnT]-i witll tlle cnem i ca i potential 
fi and the temperature T. The formula valid also for 
the Hamiltonian (|2"5|) is discussed in Appendix [B] By 
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FIG. 3: Band structures of multilayer graphenes, TV = 
2, 3, 4, 5 with 73/70 = 0.1, around the K point (taken as ori- 
gin) along the k x axis. Right panel shows a zoom-out of the 
left. Numbers assigned to curves indicate m. 



integration by parts in Eq. (|30p , we have 



(32) 



showing that the susceptibility at non-zero temperature 
is written in terms of that at zero temperature. The 
integration of x over ^ is independent of T. 

We include the impurity scattering effects by introduc- 
ing a self-energy — ir in the Green's function, i.e., iO in 
([50)) is replaced by iT. Here we simply assume the scat- 
tering rate r = H/2t to be independent of energy. 

Using the decomposition of the Hamiltonian, the mag- 
netization of the TV-layered graphene can be written as a 
summation over each sub-Hamiltonian. The contribution 
from m = is exactly equivalent to the susceptibility of 
a monolayer graphene which becomes at zero tem- 
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(33) 



Thus, the odd-layer graphene always has a large diamag- 
netic peak at zero energy. The delta-function dependence 
of Xmono agrees with the general property of the suscep- 
tibility in systems described by the fc-linear Hamiltonian, 
as discussed in Sec. IIVI 

In the presence of disorder, the delta-function is broad- 
ened into a Lorentzian with width T and the same area^i 
i.e. 



Xx 



ffvSs e 2 7 2 



67T ft 2 Tr^+r 2 ) 



(34) 



within the present model assuming a constant T. The 
shape of the peak itself depends on the model disorder 
and we may have some different manner of broadening 
in a more realistic treatment. In fact, in the mono- 
layer graphene it was shown in a self-consistent Born 
approximation 8 that x nas a much sharper peak at e = 
than the Lorentzian and also a large tail proportional to 
|e| — 1 for e 7^ 0M- In multi-layer cases effects of disorder 
are more complicated because of the presence of other 
bands. This problem is out of the scope of this work. 

The susceptibility of a bilayer graphene described by 
the Hamiltonian (fl?3")) was analytically calculated for the 
case of 73 = 0. — The expression for T = and T = is 
given by 
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with A = \N,mi where 9(t) is a step function defined by 



0(t) = 



(t > 0); 
(t < 0). 



(36) 



The susceptibility diverges logarithmically toward ep = 
0, becomes slightly positive for |£p| <J A71, and vanishes 
for \ef\ > A71 where the higher subband enters. In the 
presence of disorder, the logarithmic peak is broadened 
approximately as cx In \J e 2 F + T 2 . 

The integration of \ m Eq. (|35|) over the Fermi en- 
ergy becomes — (g v g s /2>ir){e 2 ^ 2 /h 2 ) independent of 71, 
which is exactly twice as large as that of the monolayer 
graphene ([55)1 . This arises due to the fact that the in- 
tegral of x over the Fermi energy is determined only by 
terms of the Hamiltonian matrix, proportional to k x or 
k y , and is independent of terms independent of k x and 
k y . A proof of this important property is presented in 
Sec. EH 

If we include the extra band parameter 73, the low- 
energy structure of the susceptibility (|35|) drastically 
changes due to the fine structure around the band touch- 
ing point. To demonstrate this, we numerically calculate 
X for the Hamiltonian (|23p in the case of the maximum 
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FIG. 4: Susceptibility of the sub-Hamiltonian fl 23 p in the 
multilayer graphene with Ajv jm = 2, 73/70 = 0.1 and several 
disorder strengths F. Energy is scaled in units of £ tr i g = 
O.O271. Dashed curves show plots for 73 = 0. Inset at the top 



is a zoom out of the top panel (F 
energy 71. 



0.02e t ri g ) with units of 



trigonal warping, Ajv,m = 2. Figure [4] shows the suscep- 
tibility as a function of Sf with several values of T. We 
take 73/71 = 0.1, where the Fermi line splitting occurs in 
lower than e tr i g = O.O271. For reference we also plot the 
result without the trigonal warping, 73 = 0, as a dashed 
curve. 

When we go down from high energy in the top panel 
(the smallest T), the susceptibility gradually deviates 
downward from the logarithmic dependence of 73 = 0, 
and takes a sharp dip at e = £trig- Remarkably we have 
a strong peak centered on e = 0, which is regarded as 
the effect of the linear dispersions around zero energy. 
The integral of x over the Fermi energy is almost con- 
stant — (g v g s /3Tr)(e 2 '-f 2 /h 2 ) as discussed in Sec.|lVl show- 
ing that the reduction in higher energies compensates the 
zero energy peak. As T becomes larger, the peak begins 
to cancel with the reduction in high energy and the effect 
of 73 eventually disappears when T 3> etrig- In contrast, 
the peak associated with the monolayer band m = be- 
comes broad in T but never vanishes, as shown in Eq. 
©3). 

Figure O shows x( £ f) of graphenes with layer num- 
ber from N = 1 to 5 with several disorder strengths T. 
For N > 3, insets show the contributions from each of 
bilayer-type bands. The result of odd N always contains 
a monolayer- like component, which is exactly the same 



as N — 1 and thus omitted in the inset. We can see that 
odd-layered graphenes exhibit a particularly large peak, 
which mainly comes from the monolayer-type band. A 
bilayer-like component contains a central peak due to the 
trigonal warping and a logarithmic tail in high energies, 
in accordance with Fig. Q] 

The layer-number dependence of the susceptibility in 
multilayer graphene has been studied for the graphite 
intercalation compounds. 29 This system can be viewed 
as independent multilayer graphenes bound by the in- 
tercalant layers, but the intercalants give a strong elec- 
trostatic potential along the stacking direction, leading 
to the charge redistribution among different layers4£ As 
a result the band structure and the magnetization are 
considerably different from our system with a uniform 
electrostatic potential in the vertical direction. 

In isolated multilayer graphenes realized in recent ex- 
periments, we may have some potential difference among 
layers depending on the experimental environment, and 
this can also be tuned by the external electric field as 
mentioned. In Sec. IIV1 we will show that, as long as the 
potential is not too strong to alter the entire band struc- 
ture, this does not change the qualitative feature of the 
magnetization. 



IV. DISCUSSION 

The zero energy peak in the bilayer-type subband orig- 
inates in Dirac-like dispersions appearing around four 
Fermi points. Using the known results in a bilayer^ 
we can show that the sequence of the Landau lev- 
els in the center pocket approximately becomes s = 
sgn(n)(v / 2AAr im 7'//)i/fn| with N = 0, 1, 2, • • • , and those 
in the three leg parts e = sgn(n)(\/6\N,mjyi)\/\n\, 
where I = \J h/(eB) is the magnetic length. Since 
the susceptibility is determined solely by Landau level 
energies, we compare this to the monolayer's sequence 
e = sgn(n)(-\/27 /l)-\/\n\ and obtain x from each pocket 
by substituting 7 in Eq. (|3"3"|) . We end up with 



X = 10 [ — I \„nm..- 



(37) 



except for a constant coming from the integral over the 
lower energy states. The zero-energy peak in Fig. |4] fits 
well to the Lorentzian with width T and the area of the 
delta- function (l37|) . as long as r<e t rig. 

The factor attached to Xmono becomes as large as 0.4 
when \N,m = 2 and 73/70 = 0.1, and therefore the singu- 
larity is not too small compared with that of the mono- 
layer. For A-layered graphene, a simple relation 
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leads to the summation of (f3"T|) over all the bilayer-type 
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FIG. 5: Susceptibility of multilayer graphenes with layer numbers N = 1 to 5, plotted against the Fermi energy. Results shown 
for several disorder strengths specified by constant scattering rate F. 



subbands, 



which is satisfied solely by 



X = W(N-1) 



0.1x(iY-l)Xr 



(39) 



In Fig.[5]the peak height becomes a little larger than this 
estimation due to mixing with the logarithmic tail. 

The delta-function dependence of \ in monolayer 
graphene is a characteristic property common to gen- 
eral fc-linear Hamiltonian. This can be shown using the 
scaling argument. We consider a Hamiltonian TL which 
contains only terms linear in k x and k y . We change the 
energy and wave number scales by an arbitrary factor a 
as 



e = ae, ki = aki, 



(40) 



then the Hamiltonian becomes formally identical un- 
der this transformation, since the coefficients of fc-lincar 
terms in the Hamiltonian remain unchanged. 

Going back to the definition of \ in (|30| and (|3ip. F(z) 
is scaled as 



F{z) = \P(S). 



(41) 



The function F should depend only on the coefficients 
of fc-linear terms and natural constants, and thus is in- 
variant under the scale transformation, namely we have 
F = F. With (l4"Tj) . we come up with a equation, 



F{z) 



1 



-F 



(42) 



m - ± 



(43) 



A constant A is related to the integral of the suscep- 
tibility x( £ f ) over the Fermi energy ef- From (|30|) . we 
generally have 

/oo />oo -t /> 

x(e)cfe = — Im / deeF(e + iO) = — d> dzzF(z), 
-oo J~oc 2i J c 

. (44) 

where the integral path C is a circle with an infinite ra- 
dius with anti-clockwise direction. In the present system, 
(|43f immediately gives the integral as ttA. This is an in- 
tegral of the real function x(e) and thus is real. Substi- 
tuting Eq. |43|) with real A in Eq. ([30]) , we finally obtain 
the explicit form of the zero-temperature susceptibility 
as 



X(e F ) = -Im 



.4 



= irAS(£ F )- 



(45) 



As discussed in Sec. HH the band structure in more 
realistic models has an energy gap around zero energy 
due to extra band parameters neglected in the present 
model. It was also mentioned that the external electric 
field along stacking direction opens an energy gap. One 
might think that the gap would strongly reduce the large 
diamagnetism at the band touching point. However, we 
can show within the effective mass approximation that 
the integral of susceptibility over ep is independent of 
any kind of matrix elements without k x and k y , which 
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are responsible for gap opening. This is obvious from 
the general expression (|4"4")l ; even if the Hamiltonian con- 
tains ^-independent terms in addition to /c-linear terms, 
they can be safely neglected in the integral as they are 
infinitesimal compared to \z\ on the path C. In the 
effective-mass model of the multilayer graphene we im- 
mediately conclude that the integral is independent of 
7i j l2^ 75, 76 j an d any other parameters independent of 
the wave vector. Thus we expect that the large diamag- 
netic peak is still visible even when a gap opens, while it 
may get broadened in energy by the gap width. The dia- 
magnetism in narrow gap systems is known in bismuth 45 
and recently studied for the gapped Dirac fermion4£ Any 
further discussion requires a direct computation of the 
magnetization including extra parameters, but we leave 
this for the future study. 

The integral of the susceptibility can be calculated by 
the Hamiltonian with the fc-independent terms dropped, 
and thus depends only on the band parameters associated 
with k- linear terms. In our model (|6|), the value is mainly 
determined by the dominant parameter 70, while 73 gives 
a correction at most of the order of (73/70) 2 ~ 0.01. The 
correction must be the second order in 73 because we 
can change 73 to —73 in the Hamiltonian with a uni- 
tary transformation multiplying the base on layer j by 
(— iy . As a result, the integral of x( £ f) for the bilayer- 
type Hamiltonian becomes almost twice as large as the 
monolayer's, and the summation over all the subsystems 
in A^-layered graphene becomes approximately N times 
as large as the monolayer's. 

It is instructive to derive the susceptibility start- 
ing from the Landau- level energies. In the monolayer 
graphene, the thermodynamic function f2 is given by 



with respect to h, we immediately have 



fi = -k B Tg v g s ^-p^2g(e n )ip(e n ), 



cp(s) =ln{l+exp[/3(M-£)]}, 



(46) 
(47) 



where (3 = l/fc^T, e n — sgn(n)7kj b V \n\ with hws = 
\/2j/l, and 17(e) a cutoff function which gradually decays 
to zero for \e\ > e c with cutoff energy e c . We can rewrite 
this as 

00 . 

n = - k BTg v g s ^ J2 (l-^ n0 )H(nh), (48) 



71=0 



where 

H(x) = g{y/x) In [f + 2 exp(/3/x) cosh(/3v/i)+exp(2/3/x)] , 

(49) 

with h = (Hujb) 2 - 

Expanding the integral 



-H(Q)+J2H(x+hj) 



3=1 

1 r 1 

H(x)dx-—h 2 H'(0) + -H'(oo) 



(51) 



12 L w 2 

up to the second order in h or in B. Then, we have 

= n Q + An, (52) 

where fio is the thermodynamic function in the absence 
of a magnetic field and 



An 



1 g v g s {huj B ) 2 /3exp(/3Q 
12 2nl 2 

9v9sl 2 
I2ttZ 4 



[l + exp(/3C)] 2 

(-^)*M*. 



(53) 



Applying the relation Afl = x_B 2 /2, we obtain Xmono 
given by (|3"3")) at zero temperature. 

We should note that the thermodynamic function in 
the absence of a magnetic field is given by 

1 f 1/2 

Q = -k B Tg v g s ——2] g(e n+t )<p(e n+t )dt. (54) 
27Tl n J-V* 

For contributions of states with |n| 3> 1, we can expand 
the above with respect to t and have to the lowest order 
in the field strength B 



An = 



2ttZ 2 96 
g v g s (fkv B ) 2 



(fi^) 4 le- 3 f(e n )-e- 2 f(e n )} 



lim 



2irl 2 48 s^+o 



ffvgs 1 

2tt/ 2 24 



(^° [ e - 2 /(e)- e - 1 /'( £ )]d £ 
[e- 2 f{e)-E- l f{e)]de) 
)5(e)de. (55) 



dm 

de 



This gives a "paramagnetic" susceptibility. For n = 0, 
on the other hand, the change in the thermodynamic 
potential is calculated as 



An 



2tt/ 2 8 



{hWB? 



de 



5{e)de. (56) 



The sum of these two contribution is the same as Eq. 
(f53")) . as is expected. 

In the bilayer graphene, the Landau level with 73 = 
in the region |e| <§; can be calculated from the Hamil- 
tonian as 1 ^ 



poo ph/2 00 ph/2 

I H{x)dx= I H{x)dx + Y J / H(x+hj)d. 

JO JO - =1 J-h/2 



(50) 



shui c y/n(n +1), 



(57) 



where lo c = eB/m* with m* defined in Eq. (|27|) . s = ±1, 
and n = 0, 1, 2, • • • . We have doubly degenerate levels at 
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zero energy (n = 0, s = ±1), while the spacing gradually 
becomes constant as n goes higher. In a similar but more 
complicated manner, the susceptibility is calculated as 



/ eh 



V 2m* 



2Trh 2 



(58) 

which correctly describes the logarithmic divergence 
around zero energy in the rigorous expression (|35p . as 
expected. A constant term independent of energy is miss- 
ing in Eq. (|58|) since this depends on all the low-energy 
bands which are neglected in this calculation. 

We can understand the logarithmic dependence intu- 
itively by looking into the Landau-level sequence. The 
Landau level energy can be expanded for large n as 



shw, 













(n-\ 






4) 


H 













(59) 



where the first term gives the constant interval, and the 
second gives a shift toward zero energy, which is rewritten 
as — (huj c ) 2 /(8e sn ). For sf < 0, for example, the change 
in the total energy due to the energy shift is calculated 
as 



AE = 



2Trh 2 



8e 



de 



B 2 



ffv g s 



e 2 7 2 



2 4tt ft 2 A7i 



In- 



(60) 

giving the ln|e^ | dependence of the susceptibility. 

For the Hamiltonian (|25j) containing terms propor- 
tional to fc±, the susceptibility formula ([50)) with (pH]) is 
no longer valid, since this was originally derived for sys- 
tems in which x commutes with H y and H yv 
The modified formula should be 



d 2 H/dk 2 r 



F(z) 



4ttL 2 h 2 



^ ^ tr yGTL x GTLy GTL X GTL y 

k 

—2GT~l x GH X GUy Gl~Ly 



— -GHyGH.xxG'Hy — -GHxGHyyGHx^ . (61) 

This is derived in Appendix [B] A scaling argument sim- 
ilar to the case of the monolayer graphene then gives 
F(z) ocl/z. This again leads to the logarithmic depen- 
dence of x on the Fermi energy, which coincides with (|58j) 
apart from a constant. 

The experimental measurements of the magnetization 
of two-dimensional electron systems were performed on 
the semiconductor heterostructures, by using the super- 
conducting quantum interference device (SQUID) 4 -^ or 
using the torque magnetometer i 49 i 50 ' 51 We expect that 
the detection of the graphene magnetism is also feasible 
with those techniques. 

We have studied the orbital magnetism of multilayer 
graphene with the Bernal stacking in the effective mass 
approximation. We have demonstrated that the Hamilto- 
nian and thus the susceptibility can be decomposed into 
those equivalent to the monolayer or bilayer bands. The 
monolayer-like band exists only in odd- layered graphencs 



and gives a strong diamagnetic peak at Ef — 0. The 
bilayer-like bands always exist and present a strong dia- 
magnetism in the vicinity of zero energy, unless the fine 
band structure caused by 73 is destroyed by the disorder. 
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APPENDIX A: EFFECTIVE MASS 
HAMILTONIAN 



We derive in the following the effective mass equation 
Eq. ^ describing states in the vicinity of K point in 
a multilayer graphene, by starting from the one-orbital 
tight-binding model. The following is nothing but a 
straightforward extension of the monolayer case i 13 i 36 In 
a tight-binding model, the wave function is written as 



3 L R-a, 



-^Z^BjiPBj)^ - R Bj) 



(Al) 



where j = 1, 2, • • • , N is the layer index, </>(r) is the wave 
function of the p z orbital of a carbon atom located at the 
origin, as a function of three-dimensional position r. Rx 
is the three-dimensional position of the site X, and p x 
is a two-dimensional component of Rj, parallel to the 
layer. 

In the model including hopping parameters 70, 71, and 
73 defined in Sec. HH Schrodinger's equation can be writ- 
ten as follows: For odd j, 



1=1 

3 

+73 \^B i+ x (Pa d + ri) + ^s J _ 1 (p Aj + r/)] ,(A2) 
3 

^B 3 {Pb 3 ) = -l^Y^A 3 {p B] + T l) 
1=1 

+71 [^PA J + 1 {P B] ) + VMy-i (P Bj ) ' ( A3 ) 
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For even j , 
eipAj (p Aj ) 

{p B] ) 



APPENDIX B: SUSCEPTIBILITY FORMULA 



-7o X^-Bj (^A/ ~ Tl "> 
1=1 

+71 \4>B i+ i{p Ai ) +V'Bi-i(PA 1 )] . ( A4 ) 
3 

-7o X ^ (p Bj +t-;) 



z=i 



+73 \$A j+ i (P Bj ~Tl) + IpAj-r (P Bj -Tl) . (A5) 
(=1 

Here we introduced the vectors from B site to the near- 
est neighboring A sites as T\ = a(0, 1/V3), i~2 = 
a(- 1/2,-1/2^3), and r 3 = o(l/2, -1/2^3), and we set 

"0A O = V'Bo = + 1 = 1pB N + 1 = 0. 

The states around K point can be expressed in terms 
of the slowly- varying envelope functions F A . , F Bj as 



<F A A PA i), 



(A6) 



lMP Bj ) = ^ » '' K/,; /'» (Ps ), (A7) 

where K = (27r/a)(l/3, l/v3), and CV-, are phase 
factors defined by 

C Aj =-u-\ C B] =1 (j: odd), (A8) 
C Aj = 1, C Sj = -w (j: even), (A9) 

with uj — exp(27ri/3). When x; is much smaller than the 
length scale of the envelope functions, we have 

V>x(p±"n)«e iK -^ ±T ') (l± Tl - ^)F x (p), (A10) 

with X — Aj or Bj. 



By substituting Eq. fAT]) with (|A10|) . into 
Schrodinger's equations (|A3|) and (|A5|) , we have 
for odd j, 

eF A] (p) = 7 k_F Bj (p) + 7 'k+ [F Bj _ x (p) + F B]+1 (p)] 
eF B] (p) = 7 k+F Aj (p) + 7i [F A ._ x (p) + (p)] , 

(AH) 

and for even j, 

eF Aj (p) - >yk-F Bj (p) + 7 i [F Bj ^ (p) + F Bj+1 (p)] , 
eF Bj (p) = 7 fc+F A . (p) + 7 'fc_ [F Aj _ ± (p) + (p)] , 

(A12) 

where k± = k x ± ife^ with fc^ = i fcj, = \-§^- We also 
used an identity 

]Te- iK ^(l Tf t?) = ^-ou-\0 i 1). (A13) 

If we rewrite this set of equations into the matrix form 
for a vector (F Al , F Bl , F Az , F B2 ,•••), we finally obtain 
the Hamiltonian matrix ([6]). 

The effective Hamiltonian for another valley K' = 
(27r/a)(2/3, 0) can be derived in a parallel way, while u> 
and w _1 are exchanged in Eq. (|A9|) . 



The susceptibility formula (|30|) with (|3 1 [) has been de- 
rived for the Luttinger-Kohn representation of the Bloch 
function^ and therefore in systems described by the 
Hamiltonian consisting of the free electron kinetic energy 
fi, 2 k /2m and terms linear in We derive here the sus- 
ceptibility formula which is valid in the general Hamilto- 
nian 7i(k) which includes fc-square terms in off-diagonal 
matrix elements as well as fc-linear terms. We shall con- 
fine ourselves to the case without electron-electron inter- 
action for simplicity. 

Consider the system described by the Schrodinger 
equation 



H(k+-A(r), r)^„(r) 



(Bl) 



with k = — iV and A being the vector potential. The 
thermodynamic function £1 is given by 

fl = -fcBT. 9s ^^ln{l + exp[/3( M -£ Q )]} 



1 



-/crT q s — / de ( ^ImTr- 

ys V J V ttJ e-H + iO 

x In (l + cxp[/3(/i-e)]}, 



(B2) 



where g s is the spin degeneracy and V is the system vol- 
ume. 

We consider an isotropic system and assume the vector 
potential A=(0, A), with 



B 

2iq' 



A{x) = —(e lqx -e' lqx ) 



B(x) = Bcos(qx), (B3) 



where we are going to take the long wavelength limit q — > 
0, for which the field causes the response the same as that 
due to a spatially uniform magnetic field. In the presence 
of this vector potential, the Hamiltonian changes from 
W(k) to H(k + Ak), with Ak = (0, Afc), where Ak = 
(e/h)A(x). The Hamiltonian can be expanded as 

H{k+Ak) =Ho + ^{AkHy+HyAk) + ^(Ak) 2 H yv , (B4) 

where Ho = H(k), H y = dHo/dk y , and H yy = 
d 2 H.o/dky. Note that in general fc-square Hamiltonian 
Afc does not commute with H y but does with TL y y 

Expanding the Hamiltonian up to the second order in 
the strength of the magnetic field B, we have 



Tr 



1 



e-H 



1 



d 



l 

~H a ) ~ (2iql 2 ) 2 de 
1 



Tr 



e-n Hyy 



H y .(B5) 



where I = \J h/eB is the magnetic length, 7i q = 7i(k+q), 
q = (q,Q), and we assumed that the system is transla- 
tional invariant (after the configuration average in the 
presence of impurities) . 
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We then expand 7i q up to the second order in q, to 
have 



Tr 



1 



1 



H e-H 

Y^jj-Q^^^GTi. x GTi.yGTi. x GTi. y — 2GTi. x GTi. x GTi. y GTi. y 
— -GHyGTi. xx GH y — —G'H x G'Hy y G'H x ^j ,(B6) 



with G = (e — Ho) -1 - This immediately gives the change 
of the thermodynamic potential Afl = £l(B) — O(0) with 
Eq. (|B2[) . The susceptibility \ 1S obtained by a relation 
An = -(l/2) x (B{x) 2 } = (l/4)xB 2 as 



xlmTr) GTL x G7i.yG7i. x GTLy — 2Gl~t x Gl~t x Gl~CyGl~Cy 



1 



1 



—GH y G1~l xx G1~ly — —G7i x G7iyyG7i x 



This gives Eq. (|30| with ([61]) for the multilayer graphene. 
When Afc commutes with TL V , we can simplify the for- 



-H, 



mula by noting in Eq. (|B5[) that 



Tr 



-Hi, 



Tr 



= Tr 



L-Ha 



-H % 



1 



n, 



E — Hq 

1 

e-H- 



-n. 



The susceptibility becomes 



*-&*/*'<«>(-; 



(B8) 



1 1 uTr ( GH x GH y GH x GHy} , 



(B7) giving Eq. ^3DJ> with Eq. (JH 



(B9) 
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